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I.  INTRODUCTION 


In  the  field  of  computational  mechanics  it  is  well  known  that  the  accuracy 
of  finite-difference  approximations  depends  on  the  fineness  of  the  grid 
system;  moreover,  the  accuracy  of  solutions  also  depends  upon  the  resolution 
of  solution  gradients.  For  use  of  a  uniform  grid  system  the  local  truncation 
error  in  the  difference  approximation  to  derivatives  can  be  the  largest  at 
the  location  where  the  solution  gradient  is  the  largest.  Hence,  one  can  expect 
to  obtain  greater  overall  accuracy  for  the  solution  if  more  grid  points  are  con¬ 
centrated  at  places  where  the  solution  gradient  is  very  large. 

For  a  transonic  viscous  flow  past  a  slender  projectile  with  sting,  the 
variation  of  the  solution  characteristics  in  the  direction  normal  to  the  body 
surface  is  somewhat  predictable  qualitatively  for  most  of  the  flow  region; 
consequently,  a  rather  well-suited  grid  size  distribution  in  the  normal  direc¬ 
tion  can  be  predetermined  and  fixed.  However,  the  solution  characteristics  of 
a  transonic  flow  problem  in  the  streamwise  direction  are  more  complex;  the 
location  of  nearly  normal  shocks  changes  with  the  value  of  free-stream  Mach 
number  as  well  as  with  the  number  of  iterations  used  in  the  process  of  conver¬ 
gence.  It  is  known  that  in  association  with  the  formation  of  shock  waves 
there  exist  extremely  steep  solution  gradients.  Therefore,  a  proper  distribu¬ 
tion  of  the  boundary  grid  points  can  be  very  crucial  to  the  accuracy  of  the 
computed  forces  acting  on  the  projectile  in  transonic  flows. 

For  an  axisymmetric  flow,  the  grid  system  generation  code  available  at  BRL 
can  provide  a  good  grid  network  for  the  transonic  flow  problem  if  the  boundary 
grid  points  are  properly  generated.  Presently,  the  boundary  points  required 
for  the  grid  generation  are  specified  and  clustered  in  accordance  with  the 
user's  intuition  and  experience;  the  resulting  grid  systems  have  provided 
accurate  results  for  a  number  of  different  flow  conditions.  It  seems,  how¬ 
ever,  that  the  accuracy  and  efficiency  of  the  solution  method  can  be  further 
improved  if  the  boundary  grid  points  required  for  grid  generation  are  adap¬ 
tively  determined  according  to  a  relevant  solution  gradient  distribution  along 
the  body  surface.  Moreover,  for  certain  flow  problems,  the  development  of  an 
adaptive  boundary  grid  point  generation  is  essential  for  improving  the  overall 
accuracy  of  solutions,  for  the  maximum  number  of  boundary  points  allowable  in 
numerical  simulation  is  often  limited  by  the  capacity  of  an  existing  computer 
system. 

The  objective  of  this  development  is  to  modify  the  existing  grid  genera¬ 
tion  code  GRIDGEN  so  that  the  user  will  have  an  option  to  generate  the  bound¬ 
ary  grid  points  adaptive  to  an  input  control  function.  The  resulting  grid 
generation  code  named  ADAPTGD  has  been  tested  successfully. 


II.  IMPLEMENTATION  TO  GRIDGEN 

In  modifying  GRIDGEN  for  adaptive  boundary  grid  point  generation,  efforts 
have  been  made  to  minimize  the  changes.  In  fact,  only  insertions  of  executable 
statements  are  made  to  GRIDGEN;  consequently,  the  resulting  grid  generation 
code  ADAPTGD  will  make  a  number  of  unnecessary  computations  in  the  process  of 
redistributing  the  boundary  grid  points.  A  substantial  modification  not 
related  to  the  adaptive  boundary  grid  generation  has  been  made  to  subroutine 
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RHS  ( — )  so  that  the  resulting  hyperbolic  grid  system  exhibits  more  desirable 
characteristics  for  finite-difference  computations.  The  modifications  made  to 
GRIDGEN  and  the  relevant  information  required  for  using  ADAPTGD  are  described 
and  discussed  in  the  following  sections. 

A.  Input  Data 


In  order  to  run  the  grid  generation  code  ADAPTGD,  the  user  must  provide 
two  data  cards  in  addition  to  the  regular  input  data  file  for  GRIDGEN  at  the 
very  beginning  of  the  input  data  file.  The  data  cards  provide  information  and 
instruction  for  the  redistribution  of  boundary  grid  points.  Definition  and 
comments  for  the  input  data  can  be  found  in  the  first  part  of  the  program  MAIN 
in  ADAPTGD.  A  listing  is  provided  in  Appendix  A. 

B.  Program  MAIN 

Only  insertions  of  executable  statements  have  been  made  to  the  program 
MAIN  of  GRIDGEN.  The  first  two  executable  statements  now  read  the  two  addi¬ 
tional  input  data  cards.  The  program  then  continues  as  before  to  generate 
inner  boundary  grid  points  by  calling  appropriate  subroutines.  However,  the 
generated  boundary  grid  points  will  be  overridden  if  the  user  calls  for  an 
adaptive  boundary  grid  point  generation. 

In  order  to  generate  adaptive  boundary  grid  points,  a  tape  or  file  F0R007 
which  contains  the  grid  points  to  be  redistributed  as  well  as  the  control 
function  (e.g.,  pressure)  at  grid  points  must  be  provided.  The  tape  F0R007 
must  be  written  in  the  format 

FORMAT  (1H  ,15 ,5E15.7) 
and  contains, in  order,  the  variables 

J,  XX(J),  YY(J),  XS(J),  YS  (J ) ,  FF(J) 

where  XX(J)  and  YY(J)  are  the  coordinates  of  Jth  grid  point  on  the  inner 
boundary,  while  XS(J)  and  YS(J)  are  those  of  the  outer  boundary.  FF(J)  is  the 

value  of  the  control  function  (e.g.,  computed  pressure  on  the  body)  at  Jth 

node.  The  redistribution  of  the  inner  boundary  grid  points  and  the  outer 
boundary  grid  points  are  determined  according  to  a  linear  function  of  the 
gradient  of  FF(J)  by  calling  the  subroutine  HADAPT  ( — ),  which  is  described  in 
Section  III.  Note  that  the  same  control  function  is  used  for  both  inner  and 
outer  boundary  grid  points  redistribution. 

C.  Subroutine  OUTER  ( — ) 

In  GRIDGEN  the  subroutine  OUTER  is  called  upon  to  generate  outer  boundary 
grid  points  for  elliptic  grid  generation.  If  adaptive  boundary  grid  points 
are  called  for  a  grid  generation,  then  the  call  for  subroutine  OUTER  can  be 
avoided.  In  order  to  maintain  the  versatility  of  ADAPTGD  to  reproduce  the 
capability  of  GRIDGEN  as  well  as  to  maintain  the  format  of  input  data  file 
used  in  GRIDGEN,  the  call  for  subroutine  OUTER  by  subroutine  ELPGRD  has  not 

been  modified;  however,  if  an  adaptive  boundary  grid  point  distribution  is 

called  for  a  grid  generation,  then  the  final  computation  for  XS(J)  and  YS(J) 
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in  subroutine  OUTER  is  bypassed.  Hence,  the  adaptive  outer  boundary  grid 
points  generated  earlier  are  intact  for  an  elliptic  grid  generation.  Only  one 
COMMON  statement  and  one  IF  statement  have  been  added  to  the  original 
subroutine  OUTER. 


III.  SUBROUTINE  HADAPT  (F,  X,  Y,  JI,  JE) 

The  subroutine  HADAPT  ( — )  is  developed  for  the  redistribution  of  nodal 
points  (X(J),  Y(J))  on  a  segment  of  curve  (e.g.,  inner  or  outer  boundary) 
between  nodes  JI  and  JE  according  to  the  characteristics  of  an  input  control 
function  F(J)  along  the  curve.  The  redistribution  of  the  nodal  points  can  be 
iterated  for  INO  times;  the  iteration  is  desirable  if  the  discrete  value  of 
the  control  function  is  used  for  determining  the  new  nodal  point  distribution. 
Consequently,  an  interpolation  function  subprogram  ATKN  (---)  has  been  called 
by  subroutine  HADAPT  to  provide  more  accurate  values  for  F(J)  at  the  new  loca¬ 
tion.  Note  that  INO  =  3  has  been  set  in  the  current  version  of  Subroutine 
HADAPT,  which  has  been  found  by  a  number  of  numerical  experiments  to  be  a  good 
choice. 

The  theoretical  background  for  the  subroutine  HADAPT  ( — )  is  briefly 
given  in  the  following.  Let  s  be  the  distance  measured  along  a  curve  in  xy- 
plane.  Assume  that  a  segment  of  the  curve  is  divided  into  n  elements  with 
node  at  s-j  or  (x-j,  yn-)  for  i  =  1,  2,  ...,  n+1.  Suppose  that  f(s)  is  the  con¬ 
trol  function  for  element  size  distribution  on  the  segment,  in  the  sense  that 
higher  resolutions  are  required  for  higher  gradient  of  f(s).  It  is  then 
assumed  in  Subroutine  HADAPT  ( — )  that  the  element  size  must  be  inversely 

proportional  to  a  linear  function  of  for  an  adaptive  nodal  point  distri¬ 
bution.  Therefore,  one  has 


a 


i  +  e(-^r)  • 

vds'j 


Haw. 

J 


(i) 


where  a  is  the  proportional  constant  and  3  is  the  weight  parameter  for  the 
gradient  of  control  function  f(s). 

The  proportional  constant  a  is  determined  by  equating  the  sum  of  n 
elements  to  the  total  segment  length.  Hence,  one  obtains 


a  (sn+l  "  Wj 


(2) 


With  the  input  data  f^  =  f(s-j)  and  (x-j,  y-j)  at  nodes,  the  gradient  of  f(s)  for 

a  given  element  As.  can  be  approximated  by 
J 
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(3) 


<£>•)  -  <f< 


j+i 


-  fj)/ssj>  4sj = 


-  x-i) 


1/2 


The  value  for  3  =  BW  can  be  selected  by  the  user.  It  is  clear  that  a  constant 
element  size  distribution  will  result  if  one  chooses  to  set  3  =  BW  =  0.  It  is 
mentioned  in  passing  that  the  ratio  of  the  largest  element  to  the  smallest 
element  increases  with  increasing  value  of  3.  A  value  of  3  -*■  BW  =  5.0  has 
been  set  in  the  current  version  of  the  program.  Note  that  higher  resolutions 
for  the  highest  gradient  region  can  be  achieved  by  taking  a  larger  value  of  B, 
but  then  poor  resolutions  will  result  in  the  smallest  gradient  region  if  the 
number  of  nodal  points  on  the  line  segment  of  interest  is  fixed. 

The  function  subprogram  ATKN  ( — )  called  by  Subroutine  HADAPT  ( — )  is  a 
general  purpose  interpolation  program  for  one-dimensional  problems.  Both 
subprograms  HADAPT  ( — )  and  ATKN  ( — )  which  are  placed  at  the  end  of  ADAPTGD 
package  are  also  given  in  Appendix  A  for  reference. 


IV.  OTHER  MODIFICATION  TO  GRIDGEN 

The  hyperbolic  grid  network  generated  from  GRIDGEN  for  a  projectile  has 
some  undesirable  characteristics  in  the  region  upstream  of  the  nose.  For 
example.  Figure  1  shows  rather  large  grid  sizes  upstream  of  a  projectile  nose 
and,  consequently,  poor  orthogonal  characteristics  at  the  line  of  symmetry 
result.  The  objective  of  this  modification  is  then  to  reduce  the  size  of 
those  large  grids  as  well  as  to  improve  their  orthogonal  property  at  the  line 
of  symmetry. 

In  GRIDGEN  the  subroutine  RHS  ( — )  is  called  upon  by  subroutine  HYPGRD  (- 
--)  to  provide  the  right-hand  side  of  the  reduced  finite  difference  equations 
for  hyperbolic  grid  generation.  The  subroutine  RHS  ( — )  has  been  modified, 
after  a  number  of  numerical  experiments,  by  introducing  a  weight  function  for 
controlling  the  area  of  the  grid  network.  With  the  same  input  data  of  Figure 

1,  the  modified  GRIDGEN  produces  a  hyperbolic  grid  network  shown  in  Figure 

2.  It  is  clear  that  the  modified  subroutine  RHS  ( — )  does  give  a  more 
desirable  hyperbolic  grid  network  for  projectile  aerodynamics  computations. 


V.  GRID  NETWORK  GENERATED  BY  ADAPTGD  PACKAGE 

The  modified  GRIDGEN  named  ADAPTGD  has  been  tested  successfully  for  gener¬ 
ating  a  number  of  two-dimensional  grid  networks  for  a  projectile.  In  these 
numerical  experiments  the  same  tape  file  F0R007,  which  is  given  in  Table  1, 
for  controlling  the  boundary  grid  point  distribution  has  been  employed.  The 
control  function  FF(J)  used  is  the  discrete  value  of  a  computed  pressure  coef¬ 
ficient  distribution  on  a  projectile  and  the  size  of  grid  network  chosen  is  60 
by  28.  Two  of  the  grid  networks  generated  are  presented  and  discussed  briefly 
in  the  following. 

An  elliptic  grid  network  generated  with  the  adaptive  boundary  grid  points 
is  shown  in  Figures  3(a)  and  3(b).  The  input  data  used  for  the  grid  system 
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are  given  in  Table  2.  Recall  that  only  the  first  two  lines  of  the  data  file 
are  new  to  the  data  file  for  GRIDGEN.  In  the  first  line  the  first  integer  "1" 
implies  that  adaptive  boundary  grid  points  are  called  for  the  grid  network 
generation  while  the  second  integer  "4"  indicates  that  there  are  four  segments 
of  boundary  grid  points  with  segment  end  points  fixed  to  be  redistributed. 
The  starting  point  and  the  ending  point  of  each  segments  selected  for  redis¬ 
tribution  are  specified  in  the  second  line  of  the  input  data  file.  As  indica¬ 
ted  in  the  input  data  file,  the  first  four  boundary  points  from  the  nose  of 
the  projectile  have  not  been  moved;  it  shows  that  the  user  can  choose  rather 
arbitrarily  any  segments  of  boundary  points  for  redistribution. 

The  same  input  data  have  also  been  employed  to  generate  a  hyperbolic  grid 
network  with  adaptive  boundary  points;  the  resulting  grid  network  is  given  in 
Figures  4(a)  and  4(b)  for  reference.  The  comparison  of  Figure  4(b)  against 
Figure  3(b)  shows  clearly  that  the  hyperbolic  network  has  better  overall 
orthogonal i ty  characteristics  which  can  be  advantageous  to  the  accuracy  of 
solutions  to  be  obtained  in  the  transformed  space. 


VI.  CONCLUSIONS 

The  modification  of  GRIDGEN  for  generating  grid  networks  with  adaptive 
boundary  grid  points  specified  has  been  described.  The  resulting  package  of 
computer  programs  named  ADAPTGD  has  been  tested  successfully  for  two- 
dimensional  straight  ray  grid  systems,  elliptic  grid  systems  and  hyperbolic 
grid  systems.  It  should  be  pointed  out  again  that  in  ADAPTGD  the  same  input 
control  function  has  been  employed  to  redistribute  the  inner  boundary  points 
as  well  as  the  outer  boundary  points  for  generating  an  elliptic  grid  network. 

The  modification  made  to  subroutine  RHS  ( — )  for  a  hyperbolic  grid  gener¬ 
ation  seems  to  have  made  the  hyperbolic  grid  network  competitive  to  the 
elliptic  grid  network.  An  axi symmetric  transonic  flow  past  a  projectile  is 
being  solved  with  an  elliptic  grid  system  similar  to  the  one  shown  in  Figure 
3(a)  as  well  as  with  a  hyperbolic  grid  system  similar  to  that  shown  in  Figure 
4(a).  The  results  obtained  under  the  same  conditions  and  at  iteration  number 
=  1000  show  that  both  grid  systems  give  the  pressure  coefficient  distribution 
of  the  same  order  of  accuracy.  A  straight  ray  grid  system  had  also  been  used 
for  the  flow  problem;  however,  a  number  of  difficulties  had  been  encountered, 
apparently  due  to  large  truncation  errors  resulted  from  high  skewness  of  the 
grid  system. 

For  the  application  of  ADAPTGD,  the  required  input  data  tape  F0R007  for 
adaptive  boundary  grid  point  generation  can  be  created  in  the  computer  program 
package  for  the  governing  equations  of  the  flow  problem.  For  example,  in  the 
axi symmetric  thin- layer  Navier-Stokes  code,  the  tape  F0R007  can  be  created  in 
the  main  program  of  the  package. 
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Figure  1.  A  Hyperbolic  Grid  Network  Obtained  from  Original  GRIDGEN 
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Figure  2.  A  Hyperbolic  Grid  Network  Obtained  from  Modified  GRIDGEN 
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Figure  4a.  A  Hyperbolic  Grid  Using  ADAPTGD  -  Total  Flow  Field 
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Figure  4b.  Detailed  Hyperbolic  Grid  Near  Body  Surface 
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TABLE  1.  DATA  FILE  F0R007  FOR  ADAPTIVE  BOUNDARY  GRID  POINT  DISTRIBUTION 


J 

XX  ( J) 

YY(J) 

XS  ( J ) 

YS(  J) 

1 

0. 3903E-01 

0. 5394E-05 

-0. 1800E+02 

0. OOOOE+OO 

2 

0.4903E-01 

0. 1305E-01 

-0. 1799E+02 

0. 1010E+00 

3 

0 . 9390E-01 

0. 2409E-01 

-0. 1799E+02 

0. 5755E+00 

4 

0. 1701E+00 

0. 4258E-01 

-0. 1796E+02 

0. 1381E+01 

5 

0. 2744E+00 

0. 6732E-01 

-0. 1789E+02 

0. 2485E+01 

6 

0. 4033E+00 

0.9702E-01 

-0. 1773E+02 

0. 3846E+01 

7 

0. 5532E+00 

0. 1303E+00 

-0. 1743E+02 

0 . 5414E+01 

8 

0 . 7208E+00 

0. 1661E+00 

-0. 1695E+02 

0.7132E+01 

9 

0 . 9026E+00 

0. 2032E+00 

-0. 1623E+02 

0 . 8928E+01 

10 

0. 1095E+01 

0. 2404E+00 

-0. 1521E+02 

0. 1071E+02 

11 

0. 1295E+01 

0. 2768E+00 

-0. 1390E+02 

0. 1238E+0  2 

12 

0. 1499E+01 

0. 3117E+00 

-0. 1230E+02 

0. 1385E+02 

- 

13 

0. 1703E+01 

0. 3444E+00 

-0. 1050E+02 

0. 1507E+02 

- 

14 

0. 1905E+01 

0. 3743E+00 

-0 . 8582E+01 

0. 1603E+02 

- 

15 

0. 2100E+01 

0. 4012E+00 

-0. 6629E+01 

0. 1674E+02 

- 

16 

0. 2285E+01 

0. 4248E+00 

-0. 4720E+01 

0. 1724E+02 

- 

17 

0. 2457E+01 

0. 4450E+00 

-0. 2920E+01 

0. 1757E+02 

- 

18 

0. 2612E+01 

0. 4619E+00 

-0. 1280E+01 

0. 1778E+02 

- 

19 

0. 2747E+01 

0. 4756E+00 

0. 1533E+00 

0. 1791E+02 

- 

20 

0. 2858E+01 

0. 4861E+00 

0. 1337E+01 

0. 1797E+02 

- 

21 

0. 2943E+01 

0. 4937E+00 

0. 2231E+01 

0. 1799E+02 

- 

22 

0. 2997E+01 

0. 4983E+00 

0. 2805E+01 

0. 1799E+02 

- 

23 

0. 3017E+01 

0. 5000E+00 

0 . 3016E+01 

0. 1800E+02 

- 

24 

0. 3047E+01 

0. 5000E+00 

0. 3046E+01 

0. 1800E+02 

- 

25 

0. 3141E+01 

0. 5000E+00 

0. 3142E+01 

0. 1800E+02 

- 

26 

0. 3288E+01 

0 . 5000E+00 

0. 3290E+01 

0. 1800E+02 

- 

27 

0. 3476E+01 

0. 5000E+00 

0 . 3479E+01 

0. 1800E+02 

- 

28 

0. 3692E+01 

0. 5000E+00 

0. 3697E+01 

0. 1800E+02 

- 

29 

0. 3925E+01 

0. 5000E+00 

0. 3932E+01 

0. 1800E+02 

- 

30 

0.4161E+01 

0. 5000E+00 

0. 4171E+01 

0. 1800E+02 

- 

31 

0.4391E+01 

0. 5000E+00 

0. 4402E+01 

0. 1800E+02 

- 

32 

0. 4600E+01 

0. 5000E+00 

0. 4614E+01 

0. 1800E+02 

- 

33 

0. 4778E+01 

0. 5000E+00 

0 . 4793E+01 

0. 1800E+02 

- 

34 

0. 4912E+01 

0. 5000E+00 

0. 4928E+01 

0. 1800E+02 

- 

35 

0. 4990E+01 

0. 5000E+00 

0. 5006E+01 

0. 1800E+02 

- 

36 

0. 5000E+01 

0. 5000E+00 

0. 5016E+01 

0. 1800E+02 

- 

37 

0. 5010E+01 

0. 4987E+00 

0. 5026E+01 

0. 1800E+02 

- 

38 

0. 5019E+01 

0. 4975E+00 

0. 5036E+01 

0. 1800E+02 

- 

39 

0. 5034E+01 

0. 4957E+00 

0. 5050E+01 

0. 1800E+02 

- 

40 

0. 5057E+01 

0. 4928E+00 

0. 5072E+01 

0. 1800E+02 

- 

41 

0. 5095E+01 

0 . 488  2E+00 

0. 5108E+01 

0. 1800E+02 

- 

42 

0. 5151E+01 

0 . 4814E+00 

0 . 5163E+01 

0. 1800E+02 

43 

0. 5230E+01 

0. 4716E+00 

0. 5241E+01 

0. 1800E+02 

- 

44 

0. 5338E+01 

0. 4584E+00 

0. 5346E+01 

0. 1800E+02 

- 

45 

0. 5478E+01 

0. 4412E+00 

0. 5485E+01 

0. 1800E+02 

- 

46 

0. 5655E+01 

0. 4194E+00 

0. 5660E+01 

0. 1800E+02 

- 

47 

0. 5875E+01 

0. 3925E+00 

0. 5878E+01 

0. 1800E+02 

- 

48 

0. 6141E+01 

0. 3598E+00 

0. 6143E+01 

0. 1800E+02 

- 

49 

0. 6458E+01 

0. 3208E+00 

0 . 6460E+01 

0. 1800E+02 

50 

ETC 

0. 6832E+01 

• 

0. 2749E+00 

0. 6833E+01 

0. 1800E+02 

FF  ( J) 

0 . 5133E+00 
0. 3318E+00 
0. 1522E+00 
0. 2076E+00 
0. 2075E+00 
0. 1793E+00 
0. 1396E+00 
0. 1010E+00 
0. 6669E-01 
0. 3519E-01 
0.709  IE-02 
0. 1890E-01 
0. 4226E-01 
0. 6226E-01 
0. 7971E-01 
0 . 9396E-01 
0. 1058E+00 
0. 1154E+00 
0. 1268E+00 
0. 1472E+00 
0. 1916E+00 
0. 2839E+00 
0. 3720E+00 
0. 4440E+00 
0. 4150E+00 
0. 3793E+00 
0. 3468E+00 
0. 2998E+00 
0. 2241E+00 
0. 1311E+00 
0. 5567E-01 
0. 2570E-01 
0 . 4884E-01 
0. 1325E+00 
0. 3496E+00 
0. 4080E+00 
0. 4410E+00 
0. 4734E+00 
0. 5000E+00 
0. 4621E+00 
0 . 4096E+00 
0. 3527E+00 
0. 2917E+00 
0. 2240E+00 
0. 1611E+00 
0. 1022E+00 
0.4348E-01 
0. 8539E-02 
0. 2258E-01 
0. 3760E-01 
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TABLE  2.  INPUT  DATA  FOR  GENERATING  FIGURE  3a  BY  ADAPTGD 


t — 1 

CN  H  CN 

00 

in 

VO 

r- 

00 

<T\ 

o 

1 — l 

CN 

00 

VO 

r- 

00 

av 
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1 — 1 

CN 

oo 
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H  o  o 

o 

o 

o 

o 

o 

o 

o 

•— 1 

I — 1 

I — l 

r— 1 

1 — t 

*— 1 

r-H 

r-H 

l—l 

l—l 

CN 

CN 

CN 

CN 

CU  CD  CD 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

< 

<C  o  o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

Q 

Q  o  o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

o 

<  C  o  o  o 
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o 
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o 
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o 
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o 

o 

o 

o 

o 

o 

o 

o 

o 

r- 
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o  o 


o 

VO 


00 


lD 


CN 

CN 

r- 

VO 


CN  t-i 
O  O  LD 


LD 


ro 

CN 

o 


•  I — I  i—l  LD 

O  •  CN  O  LO  • 
<Ti  O  •  •  •  iH 


LO 

r- 

*— t 

00 

1 — 1 

VO 

CN 

• 

o 

o 

o 

i — 1 

m 

• 

• 

rH 

00 

rH 

o 

r- 

• 

• 

• 

• 

• 

00 

00 

• 

o 

o 

in 

o 

I— 1 

l—l 

o 

• 

• 

• 

00 

o 

rH 

m 

o 

• 

CN 

in 

in 

r- 

m 

CN 

CN 

CN 

CN 

r- 

00 

00 

00 

rH 

00 
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vo 

o 

o 
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a> 
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00 
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• 

• 
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vo 

vo 

VO 

VO 

VO 

o 

m 

00 

in 

r-H 

rH 

• 

vo 

rH 

00 

vo 

00 

1— 1 

00 

00 

00 

o 

• 

00 

O 

CN 

00 

CN 

m 

in 
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CN 

CN 

CN 

r- 

00 

r- 

00 

00 

00 
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m 
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CN 

VO 
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LO 

o 

o 

o 

o 
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CN 

o 

o 

O 

r^* 

o 

o 

vo 

vo 

iH 
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CN 

• 

o 

O 

• 

• 

• 

VO 

• 

o 

1 
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• 

o 
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in 
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r- 

• 

00 

O 
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00 

vo 

• 

00 

o 
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o 
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o 

LO 

CN 


O  r-H 
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APPENDIX  A.  PARTIAL  PROGRAM  LISTING 


The  implementation  and  modification  made  to  the  grid  generation  package 
GRIDGEN  are  given  in  this  Appendix.  Note  that  the  additions  made  are  given 
between  the  comment  statements  CCHSU  while  the  deletion  made  in  Subroutine  RHS 
( — )  are  commented  with  CCH. 

The  Subroutine  HADAPT  ( — )  developed  for  redistribution  of  points  along  a 
curve  and  a  general  purpose  interpolation  program  ATKN  ( — )  are  also  given  in 
the  Appendix  for  references. 


4 
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nn  nonnooo  onnonoonnonnoonnn 


PROGRAM  MAIN 

COMMON  JMAX,  KMAX ,  JM,  KM,  NBOD , 
COMMON  /BOUDY/  XX (100),  YY(100), 
1  ,  T( 100 ) ,  TS ( 100 ) 

COMMON  /GRID/  X(80,60),  Y(80,60) 

CCHSU 


MAIN 

JbOD  MAIN 

XS(IOO),  YS ( 100 ) ,  SS(IUO),  S(10G)  MAIN 

MAIN 

MAIN 


COMMON/PYHSU/ IADAPT 
DIMENSION  IE( 8 ) , FF( 100) 

READ  (5,100)  IADAPT,  NSADT 

READ  (5,100)  IE( 1) ,IE( 2) ,IE( 3) , IE{ 4 ) , IE( 5) , IE( 6 ) , IE{ 7 ) , I£( 8 ) 

IADAPT=*0  FOR  REGULAR  GRID  GENERATION — GRIDGEN 1 
IADAPT=1  FOR  ADAPTIVE  BOUNDARY  GRID  POINT  GENERATION 
CONTROL  DATA  TO  BE  READ  FROM  FOR00? 

NSADT  IS  THE  NO.  OF  BOUNDARY  SEGMENTS  TO  BE  CONSIDERED  FOR  REDISTRIBUTION 
NSADT  .LE.  4  IN  THIS  PROGRAM 

IE ( K )  ARE  END  NODAL  NO.  OF  NSADT  SEGMENTS.  FOR  EXAMPLE,  THE  FIRST  SEGMENT 
IS  BOUNDED  WITH  IE(1)  AS  THE  INITIAL  NODE  AND  ENDED  WITH  IE ( 2) 

IF  NSADT  *  3,  THEN  ANY  NUMBER  CAN  BE  ASSIGNED  TO  IE (  7 )  AND  IE ( 8 ) 


MAIN 

MAIN 

MAIN 

MAIN- 

MAIN 

MAIN 

MAIN 

MAIN 

MAIN 

MAIN 

MAIN 


READ  JMAX, KMAX, IDAT,STOT, CDS  FOR  HYPGRD  SOLVER  MAIN 

IDAT  *  0  FOR  INPUT  PTS  OR  ANALYTIC  SHAPE  MAIN 

IDAT  =  1  FOR  INPUT  PTS  ON  “CARDS'*  MAIN 

IDAT  =  2  FOR  INPUT  PTS  ON  XYFILE  { TAPE7 )  MAIN 

STOT  *  -1.  FOR  CONSTANT  DS  MAIN 

MAIN 

READ  (5,80)  JMAX, KMAX, IDAT, STOT, CDS  MAIN 

WRITE  (  b , 90 )  JMAX, KMAX, I DAT, STOT, CDS  MAIN 

JM= JMAX- 1  MAIN 

KM- KMAX- 1  MAIN 

IF  (IDAT.EQ.0)  GO  TO  10  MAIN 

IF  ( IDAT. EO. 1)  CALL  INPCD  MAIN 

IF  (IDAT.EQ.2)  CALL  XYFILE  MAIN 

GO  TO  20  MAIN 

10  CONTINUE  MAIN 

MAIN 

DISTRIBUTE  POINTS  ALONG  INNER  BOUNDARY  MAIN 

WRITE  (6,110)  MAIN 

CALL  BODY  MAIN 

READ  (5,100)  NFLAG  MAIN 

IF  (NFLAG. LT.0)  GO  TO  60  MAIN 

READ  (5,100)  NCGRD , NCART  MAIN 

WRITE  (6,120)  NCGRD, NCART  MAIN 


1 3D  =  0  FOR  2D,  NO  SPIN 
1 3 D  =*  1  FOR  3D,  SYM3  HALF  SPIN 
1 3D  =  2  FOR  3D,  FULL3D  FULL  SPIN 

ISOLV  *  0  FOR  ELLIPTIC  SOLVER 
ISOLV  *  1  FOR  HYPERBOLIC  SOLVER 

READ  (5,100)  ISOLV, I3D,ND,LMAX 
WRITE  (6,70)  ISOLV, I3D,ND, LMAX 


o  o 


CCHSU 


IK  (NCGRD.GT.O)  CALL  STING  ( NCGRD) 
IF  (NCART.GT.O)  CALL  CARTB  ( NCART ) 


MAIN 

MAIN 


IF  ( IADAPT  .EQ.  0)  GO  TO  19 

11  FORMAT) 1H  ,  15 , 5E15 . 7 ) 

DO  12  1=1 ,  JMAX 

12  READ  (7,11)  J,XX( I) ,  Y  Y  ( I) ,XS( I) ,  YS  ( I) ,FF( I) 
11  =  1 

DO  13  NS=1 , NSADT 
12=11+1 
JI  =  IE(  ID 
JE=IE( 12) 

CALL  HADAPTl FF,XX,YY,JI,JE) 

IF  ( ISOLV  .EO.  1)  GO  TO  13 
CALL  HADAPT(FF,XS,YS,JI,JE) 

13  11=11+2 
19  CONTINUE 

CCHSU 

CCHSU 


JMAX=JBOD  MAIN 

JM= JMAX- 1  MAIN 

WRITE  (6,130)  JMAX  MAIN 

20  WRITE  (6,140)  MAIN 

DO  30  J= 1  , JMAX  MAIN 

WRITE  (6,150)  J ,XX( J) ,  YY( J)  MAIN 

30  CONTINUE  MAIN 

C  MAIN 

IF  ( ISOLV. EQ.0)  CALL  ELPGRD  MAIN 

IF  ( ISOLV. EQ.l)  CALL  INITIA  MAIN 

IF  ( ISOLV. EQ.l)  CALL  HYPGRD  ( STOT,CDS , IER)  MAIN 

IF  (ISOLV. EO.l. AND. IER. NE.0)  GOTO  60  MAIN 

MAIN 

FORM  3D  GRID.  LMAX  IS  CIRCUMFERENTIAL  DIRECTION  MAIN 

IF  (I3D  .EO.  2)  CALL  SPNFUL  ( I3D, ND, LMAX)  MAIN 

IF  (I3D  .EQ.  1)  CALL  GSPIN  ( I 3D, ND , LMAX)  MAIN 

60  STOP  MAIN 

C  MAIN 

70  FORMAT  (1H1,21H  ISOLV,  I3D,  ND,  LMAX, 415)  MAIN 

80  FORMAT  (3I5.2F10.0)  MAIN 

90  FORMAT  ( 1H0,27HJMAX,  KMAX ,  IDAT,  STOT,  CDS , 31 5 , 2F10 . 5 )  MAIN 

100  FORMAT  (815)  MAIN- 

110  FORMAT  ( / 1H0 , 36H++++++++++  INNER  BOUNDARY  ++++++++++)  MAIN 

120  FORMAT  (1H0.13H  NCGRD, NCART  ,215)  MAIN 

130  FORMAT  (1H0.21H  FINAL  VALUE  OF  JMAX  ,15)  MAIN 

140  FORMAT  (1H0.43H  FINAL  VALUES  OF  J,X,Y  ALONG  INNER  BOUNDARY)  MAIN 

150  FORMAT  (1H  ,I5,2F13.6)  MAIN 

160  FORMAT  (I5.6F12.5)  MAIN 

END  MAIN 
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SUBROUTINE  OUTER  ( NSEGS , IOUTD)  OUTER 

COMMON  JMAX,  KMAX ,  JM ,  KM,  NBOD,  JBOD  OUTER 

COMMON  /BOUDY/  XX(IOO),  YY(100),  XS<100),  YS(100),  SS(iOO),  S(1U0)  OUTER 

1  ,  T( 100) ,  TS ( 100)  OUTER 

COMMON  /COMP/  X(100),  Y(100)  OUTER 

COMMON  /ARRAY/  A(10Q),  B(100),  C(IOU),  D(10U),  F(1U0),  H(100)  OUTER 

CCHSU 

COMMON/PYHSU/ IADAPT 

CCHSU 

C  OUTER 

C  THIS  PROGRAM  FORMS  AN  OUTER  GRID  BOUNDARY  USING  CONTIGUOUS  CUBIC  OUTER 


C  SEGMENTS.  NUMBER  OF  SEGMENTS  IS  NSEGS.  POINT  AND  SLOPE  ARE  INPUT  AT  OUTER 


C  THE  ENDS  OF  A  SEGMENT.  SLOPE  IS  AN  ANGLE  IN  DEGRESS.  PARAMETRIC  OUTER 

C  CUBICS  USED  TO  PERMIT  ANY  SLOPE  (THETA  *  90,  -90,  ETC).  INITIAL  OUTER 

C  LOGIC  DETERMINES  CUBIC  COEFFICIENTS  OF  EACH  SEGMENT.  REMAINING  OUTER 

C  LOGIC  DISTRIBUTES  POINTS  ALONG  OUTER  BOUNDARY  USING  ARC  LENGTH  OUTER 

C  AS  DISTRIBUTION  FUNCTION.  THUS  TWO  PARAMETRIC  VARIABLES  ARE  USED.  OUTER 

C  FINDING  X.Y  SO  CUBIC  SEGMENTS  CAN  BE  DISJOINT  IN  SLOPE  IS  MESSY  OUTER 

C  A  SINGLE  SPLINE  INTERPOLATION  CANNOT  BE  USED  OVER  THE  COMBINED  OUTER 

C  SEGMENTS  BECAUSE  OF  POSSIBLE  SLOPE  DISCONTINUITY.  OUTER 

C  OUTER 

DIMENSION  JA( 8 ) ,  JB(8)  OUTER 

DIMENSION  CA0 ( 8 ) ,  CA1(8),  CA2(8),  CB0(8),  CB1(8),  Cb2(8),  CARC(8)  OUTER 
DO  110  N*l, NSEGS  OUTER 

C  POINTS  AND  SLOPES,  -90  . LE .  THETA  . LE .  90  DEGREES  USED  OUTER 

READ  (5,240)  X0 , Y0 , XI , Yl ,TH0 ,TH 1  OUTER 

WRITE  (6,250)  X0 , Y0 ,X1 , Yl ,TH0 ,TH1  OUTER 

RTH0*0 . 01745329 2*TH0  OUTER 

RTH 1=0. 01 745329 2*TH1  OUTER 


190  CONTINUE  OUTER 

CCHSU 

IF  (IADAPT  .EQ.  1)  GO  TO  220 

CCHSU 

C  OUTER 

C  FORM  PARAMETRIC  ARRAYS, FROM  DISTRIBUTED  PARAMETRIC  ARRAY,  OUTER 

C  USE  IT  TO  DETERMINE  X,Y  WITHIN  A  OUTER  SEGMENT  CURVE.  OUTER 

C  SPLINE  REQUIRES  ABOUT  5  POINTS  IN  AN  INTERVAL  OUTER 

S(1)=0.  OUTER 


END  OUTER 


f 
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SUBROUTINE  RHS  (K,SY,M,N)  RHS 

COMMON  JMAX,  KMAX ,  JM,  KM,  NBOD,  JBOD  RHS 

COMMON  /TRXE/  XXSI(IOO),  YXSI (100),  XADA(IOO),  YADA(IOO)  RHS 

COMMON  /GRID/  X(80,60),  Y(80,60)  RHS 

COMMON  /VOL/  V(100),  VSTAR(IOO),  DS(60)  RHS 

DIMENSION  SY( 160 )  RHS 

C  FILL  RIGHT-HAND-SIDE  VECTOR  FOR  HYPER  SOLVER  RHS 

MSY=*1  RHS 

DO  10  J*1 ,  JMAX  RHS 

S Y( MS Y) =XXS I(J)*X(J,K) + YXSI ( J ) *  Y ( J , K )  RHS 

S  Y( MS Y+l ) =- ( YXSI (J) ) *X( J , K ) +XXSI ( J)*Y(J,K)+V(J) +VSTAR( J )  RHS 

MSY=MSY+2  RHS 

10  CONTINUE  RHS 

C  FRACTION  OF  EXPLICIT  DISSIPATION  THAT  IS  NOT  PUT  IN  RHS 

C  IMPLICITLY  RHS 

EPS= .02  RHS 

C  CAUTION  —  EPS  IS  ALSO  SET  IN  SUBROUTINE  MATRX  RHS 

CCH  JMM=JM— 1  RHS 

CCH  DO  20  J=3 , JMM  RHS 

CCH  MS  Y=2  * J- 1  RHS 

CCH  SCALE=EPS*SQRT( XADA( J) **2+YADA( J) **2)  RHS 

CCH  XX=X( J-2,K)-2.*(X( J-1,K)+X( J+1,K) )+2.*X< J,K)+X( J+2,K)  RHS 

CCH  YY=Y(J-2,K)-2.*(Y(J-1,K)+Y(J+1,K))+2.*Y(J,K)+Y(J+2,K)  RHS 

CCH  SY ( MS Y ) =SY( MS Y ) —SCALE* ( XXS I ( J ) *XX+  YXS I  (  J )  *  Y Y)  RHS 

CCH  SY(MSY+1)=SY(MSY+1) —SCALE* ( — YXSI ( J ) *XX+XXSI ( J ) *  YY )  RHS 

CCH20  CONTINUE  RHS 

CCHSU 

L=1 

DO  30  J=1 , JMAX 
AJ=J 

SC ALE* EPS* SORT ( XADA( J) **  2+ YADA( J) **  2) 

IF  (J  .EQ.  1)  GO  TO  31 
IF  (J  .EQ.  JMAX)  GO  TO  32 
XX=X( J+1,K)-2.*X( J,K)+X( J-1,K) 

YY*Y( J+1,K)-2.*Y( J,K)+Y(  J-1,K) 

GO  TO  35 

31  XX=2.*X(J+1,K)-2.*X(J,K) 

YY=0. 

GO  TO  35 

32  XX=0. 

YY=-2.*Y( J,K)+2.*Y( J-1,K) 

35  CONTINUE 

SY(  L)=SY( L) -SCALE* ( XXSI( J ) *XX+ YXSI ( J ) *YY) 

CHEN=( 1.-60. /AJ) 

IF  (J  .LT.  4)  CHEN=( l.-( 100.-AJ*10. )/AJ) 

SY( L+l ) =SY( L+l ) -SCALE* ( - YXSI ( J ) *XX+XXS I ( J ) * Y Y ) *CHEN 
30  L=L+2 
CCHSU 


RETURN 

END 
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RHS 


SUBROUTINE  HADAPT{ FF ,XX , YY , JI , J£) 

C  FF( I )  IS  A  FUNCTION  OF  BOUNDARY  COORDINATE  SUCH  AS  PRESSURE  DISTRIBUTION 
C  ON  A  PROJECTILE;  ITS  GRADIENT  IS  USED  FOR  REDISTRIBUTION  OF  THE 

C  BOUDARY  GRID  POINTS 

C  XX(I)  AND  YY( I )  ARE  (X,Y)  OF  THE  GIVEN  BOUDARY  POINTS 
C  JI  IS  THE  INITIAL  NODAL  NUMBER  WHILE  JE  IS  THE  END  NODAL  NUMBER 
C  NOTE  THAT  THE  BOUNDARY  POINTS  BETWEEN  I»JI  AND  I»JE  ARE  TO  BE  REDISTRIBUTED 
DIMENSION  FF( 100) ,XX( 100) , YY( 100) rSX( 100) fSY( 100) ,SF< 100) f FO( 10 
10) ,DS( 100) ,TX( 100) ,TY( 100) ,TF( 100 ) , DT( 99 ) , G< 99 ) , W( 99 )  ,  XO( 100) 

N-JE-JI 

NE-N+1 

IF  (N  .GT.  99)  GO  TO  120 

C  INO  IS  THE  NUMBER  OF  ITERATION  FOR  REDISTRIBUTION 
I  NO=  3 
ITER»0 

DO  10  J=1,NE 
K=JI-1+J 
FO( J)=FF( K) 

XO(  J)=*XX(  K) 

SF(  J)=*FF(  K) 

SX(J)»XX(K) 

10  SY( J) =YY( K) 

TF ( 1)»SF( 1) 

TX  (  1)»SX( 1) 

TY( 1)«SY( 1) 

TFtNE)-SF(NE) 

TX(NE)»SX(NE> 

TY( NE) aSY( NE) 

20  CONTINUE 

C  BW  IS  THE  WEIGHT  PARAMETER  FOR  THE  GRADIENT  OF  FF 

C  THE  EFFECT  OF  THE  GRADIENT  ON  REDISTRIBUTION  INCREASES  WITH  INCREASING  BW 
BW»5.0 
ST»0. 

WT=*0. 

DO  30  J=1,N 
J1=J+1 

DS( J)=SORT< (SX( Jl )-SX( J) ) **2+(SY( J1)-SY( J) )**2) 

ST=ST+DS( J) 

G< J)»(SF( J1)-SF( J) )/DS( J) 

W( J)»l./( l.+BW*ABS(G( J) ) ) 

30  WT=WT+W(J) 

ALPHA»ST/WT 
DO  40  J=1,N 
40  DT( J) =ALPHA*W( J) 

TJ=0. 

SK=0. 

K=0 


% 


30 


f 


DO  80  J= 1 , N— 1 
J1=J+1 
TJ=TJ+DT( J) 

50  K=K+1 

SK=SK+DS( K) 

IF  (TJ  .GE.  SK)  GO  TO  50 
SK=StC-DS(  K) 

R«(TJ-SK)/DS(K) 

K1=K+1 

TX( J1)=SX( K)+R*(SX( K1)-SX( K) ) 

TY( J1 )»SY( K)+R*(SY( K1)-SY( K) ) 

TXI=TX( Jl) 

TF( Jl)=ATKN(XO,FO,NE,2,TXI) 

K=K-1 

80  CONTINUE 

81  FORMAT (  1H1  ,  ' ***  BW  =»  ’  ,  F8.4) 

WRITE( 6,81)  BW 

82  FORMAT ( 1H0,,J,.9X,'DS,,9X,,DT,»9X,,SF',9X,'TF,«9X,'SX,,9X,TX) 

83  FORMAT ( 1H  ,I3,6E12.4) 

WRITE ( 6,82) 

DO  84  J=1,N+1 

84  WRITE ( 6,83)  J , DS( J) , DT( J) ,SF( J) ,TF( J) ,SX( J) ,TX( J) 

ITER= ITER+ 1 

IF  (ITER  .GT.  INO)  GO  TO  100 
DO  90  J=2,N 
SF( J) =TF( J) 

SX( J)*TX( J) 

90  SY ( J) =TY( J) 

GO  TO  20 

100  DO  110  J»2,N 
K-JI-l+J 
XX( K) =TX( J) 

110  YY(  K) =TY( J) 

RETURN 

120  WRITE( 6,130) 

130  FORMAT ( 1H1,'****  NO.  OF  POINTS  EXCEEDS  THE  DIMENSION  ****') 
STOP 
END 
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FUNCTION  ATKN(X,Y,N,K,XI) 

AITKEN  INTERPOLATING  FUNCTION 

X(I),I=1,N, —  INDEPENDENT  VARIABLE  IN  ASCENDING  OR  DESCENDING  ORDER 
Y(  I ) — TABLE  OF  DEPENDENT  VARIABLE 

X —  DEGREE  OF  INTERPOLATION  DESIRED;  K  .LE.  12 

XI—  X-VALUE  WHERE  THE  INTERPOLATION  IS  DESIRED 
DIMENSION  X(N) , Y(N) ,XX( 13 ) , YY( 13) 

DATA  KMAX/12/ 

IF  (K  .GT.  KMAX  .OR.  K  .LE.  0)  GO  TO  300 
K 1 3  K+ 1 


IF  (X(N)-X(l))  100,10,10 
10  IF  (XI-X(l))  20,20,30 
20  LL»0 

GO  TO  200 

30  IF  (X(N)-XI)  40,40,50 
40  LL-N-Kl 
GO  TO  200 
50  LL=1 


LU3N 

60  IF  (LU-LL-1)  180,180,70 
70  LI*( LL+LU)/2 

IF  (X(LI)-XI)  80,80,90 
80  LL»LI 
GO  TO  60 
90  LU3LI 
GO  TO  60 

100  IF  (XI-X(l))  120,20,20 
120  IF  (X(N)-XI)  130,40,40 
130  LL=1 
LU»N 

140  IF  (LU-LL-1)  180,180,150 
150  LI»(LL+LU)/2 

IF  (X(LI)-XI)  160,170,170 
160  LU-LI 

GO  TO  140 
170  LL=LI 

GO  TO  140 
180  LL»LL-<  Kl  +  D/2 

IF  (LL)  20,200,190 
190  IF  ( LL+K1-N)  200,200,40 
200  DO  210  1-1, K1 
I13LL+I 

XX(  I)=*X(  I1)-XI 
210  YY ( I ) =  Y ( 1 1 ) 

DO  220  I»1,X 
DO  220  J»I,K 

220  YY( J+1)»(1.0/(XX( J+1)-XX( !))>*( YY< I)*XX< J+l) 
1-YY( J+l )*XX( I) ) 

ATKN* YY ( Kl) 

RETURN 

300  WRITE( 6,301)  K 

301  FORMATt ' 1 POLYNOMIAL  OF  DEGREE  K  »',I3,2X,'IS 
1  INCORRECT  FOR  THE  FUNCTION  SUBPROGRAM  ATKN '  ) 

RETURN 

END 
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